function y = dlsim(sys,u,v,t)
% dlsim(sys, u, v, t) returns the time response of the lti model 
% sys', where sys is of the form y(z) = B(z)/A(z) * u(z) und sys' 
% is y(z) = B(z)/A(z) + 1/A(z) * v(z), that is, a _d_isturbance v
% (hence dlsim) is added in every step. We need this specific 
% model structure because it is the model structure identified by
% the basic least squares model identification algorithm, that we
% have programmed and want to test.
% At the moment t is ignored.

[num,den]=tfdata(sys,'v');

if(abs(den(1)) < 5*eps)
    error('first entry of denominator is singular. this case is not handled yet.')
end

if(any(size(num) ~= size(den)) || numel(num) ~= length(num)) 
    error('I was expecting numerator and denominator vectors returned by tfdata to always be the same length and one-dimensional. This assumption was just violated.')
end

n=length(num);
%checks whether sizes of u and t and v match
if(not(all(size(u) == size(t)) && all(size(t) == size(v))) || numel(u) ~= length(u) )
    error('Inputs in wrong format.')
end

% make vector dimensions parallel 
if(min(size(num) .* size(u)) ~= 1)
    num = num.';
    den = den.';
end

N = length(u);

y = zeros(size(u));

%order is important here
num = num / den(1);
den = den / den(1);

for i=1:(n-1)
    y(i) = sum(y(1:i-1) .* -den(i:-1:2)) + sum(u(1:i) .* num(i:-1:1)) + v(i);
end

for i=n:N
    y(i) = sum(y((i-n+1):(i-1)) .* -den(end:-1:2)) + sum(u((i-n+1):i) .* num(end:-1:1))  + v(i);
end

end